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Abstract 

Recent work on the analysis of Iced airfoils 
and wings Is described. Ice shapes for multi- 
element airfoils and wings are computed using an 
extension of the LEWICE code that was developed 
for single airfoils. The aerodynamic properties 
of the Iced wing are determined with an Inter- 
active scheme In which the solutions of the 
Invlscld flow equations are obtained from a panel 
method and the solutions of the viscous flow 
equations are obtained from an Inverse three- 
dimensional finite-difference boundary- layer 
method. A new Interaction law Is used to couple 
the Invlscld and viscous flow solutions. 

The newly developed LEWICE multielement code 
Is applied to a high-lift configuration to 
calculate the Ice shapes on the slat and on the 
main airfoil and on a four-element airfoil. 

The application of the LEWICE wing code to the 
calculation of Ice shapes on a MS-317 swept wing 
shows good agreement with measurements. The 
Interactive boundary-layer method Is applied to a 
tapered Iced wing In order to study the effect of 
Icing on the aerodynamic properties of the wing 
at several angles of attack. 

1 .0 Introduction 

In recent years there has been considerable 
research activity In the area of aircraft Icing 
to combat the adverse effects of leading-edge Ice 
formation on fixed and rotary wing aircraft and 
on engine Intakes. Computational work and related 
experimental studies have been Initiated and are 
being carried out under the NASA Aircraft Icing 
Research Program to develop and validate a series 
of mutually compatible computer codes to predict 
the details of an aircraft Icing encounter. 7 
The papers presented each year at the the AIAA 
Aerospace Sciences Meeting and the papers pre- 

sented In this symposium show that Indeed much 
progress has been made In this area. 

In this paper we report a summary of our prog- 
ress In predicting Ice shapes on airfoils and 
wings and In determining the effect of Ice forma- 
tion on aerodynamic performance degradation. For 
airfoil flows, our research has led to Improve- 

ments In the LEWICE code 2 for predicting leading- 
edge Ice formation 3 and to the development of an 
Interactive boundary-layer (IBL) method 4 for 
determining the Increase In drag and loss of lift 

of airfoils 5 * 6 and helicopter blades 7 due to 

Icing. This capability for predicting Ice shapes 
on airfoils has also been extended by the authors 
to Include airfoils with slats, and very recently 
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to multielement airfoils so that Ice shapes on the 
main airfoil and on the flap can be computed as 
well as on the slat. This will permit the effects 
of Icing on high- lift configurations to be com- 
puted using the Interactive boundary-layer method 
recently developed by Cebecl et al. 8 

For wing flows, our research followed a similar 
path, concentrating on the development of (1) a 
three-dimensional version of the LEWICE code, (2) 
a three-dimensional Interactive boundary-layer 
(IBL) method for Iced wings, and (3) the coupling 
of the IBL method to the LEWICE code to determine 
the Ice shapes and their effects on lift, drag 
and moment coefficients for wing flows. 

The progress to date for airfoil flows Is des- 
cribed In several papers. For this reason the 
present paper concentrates mainly on three- 
dimensional flows and describes the extension to 
wing flows of the combined LEWI CE/IBL procedure 
developed for airfoils. Section 2 describes the 
method for calculating Ice shapes on the leading 
edge of a wing and presents a comparison between 
calculated and experimental results. Section 3 
describes the Interactive boundary-layer method 
for computing three-dimensional flows on Iced 
wings. In addition, this section presents the 
results obtained from the application of this 
method to a NASA MS 317 tapered wing with Ice and 
to an unswept NACA 0012 wing without Ice. In each 
case, the Invlscld and viscous flow equations are 
solved Interactively to determine the Increase In 
drag due to Ice” and to compare the calculated 
pressure coefficients with measured values. Sec- 
tion 4 presents recent results obtained for multi- 
element airfoils and Is followed by concluding 
remarks. 

2.0 Extension of LEWICE to Wings and 
Its Validation 

The extension of the LEWICE airfoil code to 
wings Is not so straightforward. There are sev- 
eral possible approaches that can be pursued. In 
each approach the flowfleld calculations should 
be performed In principle using a three-dimen- 
sional Invlscld method, and the impingement pat- 
tern of the water droplets on the surface should 
be determined by performing trajectory calcula- 
tions for the three components of the velocity 
obtained from the Invlscld method. The heart of 
the LEWICE code, however, Is the third module that 
contains the quasi-steady-state surface heat 
transfer analysis in which mass and energy equa- 
tions are solved for a two-dimensional flow In 
order to determine the Ice shape and size. The 
extension of this module to three-dimensional 
flows would require the heat balance equation, 
developed for airfoil flows, to be modified to 
wing flows. And, as discussed In Ref. 9, this 
can only be done with the help of experimental 
data that presently do not exist. As a first 
step, It Is best to leave the heat balance In Its 
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two-dimensional form and assume It to apply to a 
three-dimensional body expressed In an equivalent 
two-dimensional form. One approach, followed by 
Potapczuk and Bldwell, 1 ® Is to perform the tra- 
jectory calculations for a three-dimensional flow- 
field and apply them along the streamlines on the 
wing. Another approach, followed In Ref. 9, Is 
to approximate the 3-D wing by an equivalent yawed 
Infinite wing at each spanwlse station. In this 
case, the flowfleld Is calculated by a three- 
dimensional panel method and the particle trajec- 
tories calculated for flow normal to the leading 
edge subject to the Infinite swept-wlng assump- 
tion. Another approach Is to apply the LEWI CE 
airfoil code to the streamwlse cross-section of 
the wing. This approach has at least two alterna- 
tives, one of which Is described In this paper. 
The accuracy of these three approaches and others 
depend on the angle of attack and the spanwlse 
location of the airfoil section, and’ they require 
a careful evaluation through comparisons with 
experimental data. 

2.1 Comparison of Measured Ice Shapes and Pre - 
dictions Obtained with Yawed Wing Approximation 

The calculated results obtained with the exten- 
sion of the LEWI C E airfoil code to wing flows by 
the method of Ref. 9 are shown In Figs. 1 and 2 
together with the experimental results 11 on an 
HS-317 swept wing. A summary of test conditions 
used In the calculations are given In Table 1. 
Additional studies for other test conditions are 
In progress and will be reported later. The calc- 
ulated Ice shapes In Figs. 1 and 2 were obtained 
for the untapered wing with a MS-317 airfoil sec- 
tion defined streamwlse with a sweep angle of 30° 
and an aspect ratio of six. All trajectory and 
Ice accretion calculations were carried out for 
Invlscld flow computed on the mld-semlspan section 
where the spanwlse pressure gradient was negligi- 
ble. All calculations were performed for one time 
step to save computer time, which Is approximately 
7 minutes per run on the Cray computer. The In- 
crease In time, In comparison with the two- 
dimensional case, Is primarily due to the trajec- 
tory calculations where, despite the yawed Infi- 
nite wing approximations, the computation of the 
off-body velocities Involves repeated large matrix 
multiplications In which all wing panels are 
represented . 

Figure la shows a comparison of measured and 
calculated Ice shapes for Run 8, which corresponds 
to T» » 0°, a = 2°, t * 390 sec. As can be seen, 
the agreement between measured and calculated 
results Is remarkably good. The calculated 
results for a calculation time of 1164 sec and for 
Too » 0° F a nd a = 2° (Run 11) are shown In Fig. lb 
and Indicate reasonable agreement with measure- 
ments despite the one time step used In the calc- 
ulations. It Is expected that the Ice growth will 
have some effect on the velocity field and on the 
calculated droplet Impingement. A comparison of 
predicted and measured Ice shapes obtained for 
Tgo = 0°F at a * 8 6 for t = 390 and 1164 sec. (Runs 
9 and 10) are shown In Figs, lc and Id, respec- 
tively. The agreement Is again reasonable, keep- 
ing In mind that only one time step was used In 
the calculations. 

The next set of studies was conducted for a 
slightly higher freestream temperature of T® 
* 1 5°F , representing an Icing condition for which 



Fig. 1. Comparison of calculated (solid lines) 
a n d meas u red (da shed 1 In es) 1c e shapes. HI me 1 c e : 
(a) Run 8, (b) Run 11, jc) Run 9, (d) Run 10. 

a mixed Ice growth was observed. Run 7 In Fig. 2a 
for a = 2“ and t = 390 sec". Indicates good agree- 
ment between experiment and theory, except for 
some deviation on the upper surface. The results 
In Fig. 2b at the large time step of t = 1164 sec. 
(Run 1) are more or less In agreement In predict- 
ing the amount of Ice accumulated, but they differ 
In predicting Its shape. It Is known from two- 
dimensional calculations that a large nu mbe r of 
relatively short time steps are needed to predict 
horn-shaped Ice for glaze Ice. Since the mixed 
Ice formation tends toward glaze Ice shapes for 
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Mg. 2. Comparison of calculated (solid lines) 
and measured (dashed lines) Ice shapes. Mixed 
Ice: (a) Run 7, (b) Run 1 , (c) Run 6. 


Table 1 . 

Test Conditions for 

MS-317 I< 

:e Accretli 

Exper J 

Iment of 

Ref. 0, Vo, * 

1 50 mph , 

d = 20 pm 



LWC =1.03 gm- J . 

Run 

iii 

t 

I sec ) 

(deg) 

<k s /c)l 

1 

15 

1164 

2.0 

0.00192 

6 

15 

1164 

8.0 

0.00192 

7 

15 

390 

2.0 

0.00192 

8 

0 

390 

2.0 

0.00127 

9 

0 

390 

8.0 

0.00127 

10 

0 

1164 

8.0 

0.00127 

11 

0 

1164 

2.0 

0.00127 


large times, It Is not surprising that one time 
step calculation Is not sufficient to predict the 
actual growth of the Ice shapes. Similar consents 
apply to Fig. 2c, where comparisons are for a 
large time step of t - 1164 sec. (Run 6), but at 

a = 8°. 

2 • 2 Comparison of Measured Ice Shapes and Pre - 
dictions Obtai ned with Strip Theory Approximation 

Additional calculations were also performed 
with the LEW1CE code to determine the Ice shapes 
on the leading edge of the MS-317 swept wing dis- 
cussed In the previous subsection. This time we 


used the strip theory approximation rather than 
the yawed wing approximation. Me calculated the 
three-dimensional velocity field from the panel 
method and used the velocity distribution In the 
LEWI CE code for the streanwlse airfoil section. 
Figure 3 shows a comparison between the Ice shapes 
computed with strip theory (2-D) and yawed wing 
(3-D) approximations together with the measured 
Ice shape for run 11. As can be seen, both calc- 
ulated Ice shapes, at least for this run, agree 
reasonably well with experimental data. Addi- 
tional studies are underway to further Investigate 
the differences between the two procedures. 



Fig. 3. Comparison of calculated Ice shapes with 
experimental data. 2-D represents the Ice shape 
with strip theory approximation and 3-D that with 
yawed wing approximation. 

2 * 3 The Role o f Mind Tunnel Effect on the Calcu - 
lation of Ice Shapes 

In general, Infinite yawed wing conditions 
apply to the mld-semlspan section of wings with 
an aspect ratio greater than about five. This 
approximation becomes progressively less accurate 
as the tip or the root of the wing Is approached, 
but In most Instances It can still provide reason- 
able answers. A point to remember about the use 
of this approximation with finite aspect ratio 
wings Is that although the flow may have the 
desired characteristics, Its lift is always less 
than the lift of a wing with Infinite aspect 
ratio. This may lead to problems In comparing 
calculations with experimental data unless the 
aspect ratio or the pressure distribution Is also 
given. If the pressure distribution Is not avail- 
able, the given angle of attack may not properly 
represent the experimental conditions. Similar 
problems may also arise In simulating wind-tunnel 
conditions by calculating the corrected Incidence 
and lift coefficient in free air, because the 
trajectories in the two cases may be far from 
Identical. One solution to the wind-tunnel prob- 
lem, which may be the only acceptable solution for 
a swept wing spanning the tunnel, Is to calculate 
the flowfleld about the wing In the presence of 
the tunnel walls. 

The comparisons between the calculated and 
experimental Ice shapes presented In Subsections 
2.1 and 2.2 were obtained for the Icing conditions 
and angle of attack given In Table 1. Care, how- 
ever, Is required to perform the calculations as 
closely as possible to the stated experimental 
conditions. Even though the atmospheric Icing 
conditions are properly defined In the LEWICE 
calculations, the angle observed In the wind 
tunnel together with the wing aspect ratio may 
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need to be different when the flowfleld calcula- 
tions are performed with the panel method for a 
free air model . 

To Investigate this possibility further, we 
have calculated the pressure distributions for a 
constant chord wing with MS-317 streamwlse sec- 
tions and with a 30° sweep, having a finite aspect 
ratio In free air and spanning the side walls of 
a wind tunnel, using the panel method of Ref, 12. 
The free air model was chosen to have an aspect 
ratio of 6 In order to reduce the root and tip 
effects at the mld-semlspan location. It was 
found that the angle of attack of the free-air 
model had to be Increased to 4° for the pressure 
distributions to match the experimental data of 
Bldwell* measured at 2° angle or attack In the 
wind tunnel. The 8° angle-of-attack case In Table 
1 required an Increase of 3.5° angle of attack In 
free air to obtain satisfactory agreement with the 
pressure distributions. Since there Is some doubt 
about the flowflelds being matched at the widely 
differing angles of attack In the wind tunnel and 
In free air, the flow was calculated about the 
wing In the wind tunnel. This requires additional 
paneling of the tunnel floor, celling and one 
sidewall, while taking advantage of one plane of 
symmetry. Figure 4 shows the calculated and mea- 
sured pressure distributions at a = 2° In the 
wind tunnel compared with results from the calcu- 
lations In free air at a = 2 6 and 4°. Agreement 
between the experimental data and the calculations 
for the wing In the wind tunnel Is very good, 
considering that the calculated pressure distribu- 
tion corresponds to Invlscld pressure distribution 
and does not Include any viscous effects. As can 
be seen, matching of free-air calculations with 



Fig. 4. Computed and experimental pressure dis- 
tributions for the MS-317 airfoil. — denotes 
results for a * 2°, free air, — a = 4°, free 
air. Symbol o denotes experimental results for 
a = 2°, wing tunnel and o denotes calculated 
results for o = 2°, wind tunnel. 


wind-tunnel data Is a trial and error process. 
Studies are underway to extend these calculations 
to Include a wing with Ice. This Is relatively 
easy, except for the longer computing times for 
the particle trajectories resulting from the large 
number of panels used In the calculation of the 
Invlscld flowfleld, 

3.0 Three-Dimensional Interactive Boundary- 
Layer Method 

A complete description of the three-dimensional 
Interactive boundary-layer method Is described In 
Ref. 9 and Is presented here for completeness. 

As In two-dimensional flows, the Interactive 
method for three-dimensional flows Is based on 
the solutions of the Invlscld and boundary-layer 
equations. An Interface program, Illustrated by 
Fig. 5, Is placed between the Invlscld and three- 
dimensional Inverse boundary-layer methods to pro- 
cess the geometry and Invlscld velocity data for 
Input to the boundary-layer program. The basic 
Input to this program Is (1) the definition of the 
wing configuration that Is used by a geometry sub- 
routine to construct a nonorthogonal coordinate 
system and (2) the associated geometrical param- 
eters, which consist of the geodesic curvatures 
and metric coefficients needed In the boundary- 
layer calculations. Some of the generated data 
Is used later In a velocity subroutine to deter- 
mine the Invlscld velocity components at the 
boundary-layer grid points and to transform the 
Invlscld velocity components on the surface, calc- 
ulated In a Cartesian coordinate system. Into the 
boundary-layer coordinate system. This operation 
consists of calculating dot products of velocity 
vectors as well as chordwlse and spanwlse Interpo- 
lation. Further velocity and geometry data pro- 
cessing Is carried out In a subroutine that sepa- 
rates the generated Information Into upper and 
lower surfaces of the wing for boundary- layer 
calculations . 



Fig. 5. The Interactive boundary-layer method. 

The above procedure Is appropriate to wings 
without Ice and has been used to compute transonic 
flows on wing/body configurations' 8 where, since 
the wing leading edge was free of Ice, there was 
no difficulty In generating solutions near the 
attachment line by constructing the nonorthogonal 
coordinate system and computing the geometrical 
properties of the wing. For a wing with Ice, gen- 
eration of the boundary-layer solutions near the 
leading edge can pose problems since the geodesic 
curvatures and metric coefficients must be deter- 
mined for an Irregular surface. In addition, the 
formulation of the Interactive boundary- layer 
method developed for Iced airfoils must take 
account of the three-dimensional nature of the 
flow. Thus, It Is necessary to make changes In 
the strategy for solving the three-dimensional 
boundary-layer equations for an Iced wing. These 
are considered below. 
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3.1 Boundary-Layer Equations 


The three-dimensional boundary-layer equations 
for a nonorthogonal coordinate system are given 
In several references. With Reynolds stresses 
modeled by the eddy-viscosity concept, they can 
be written as, 


~ (uh^ sine) ♦ ^ (wh ] sine) ♦ (vh 1 h 2 sine) * 0 
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Here x denotes the coordinate along the lines 
formed by the Intersection of the wing surface and 
planes representing constant percent semispan; z 
Is the coordinate along the constant percent 
chordllnes that generate the wing surface, with 
chord defined as the maximum length line between 
leading edge and trailing edge. The third coord- 
inate y denotes the direction normal to the wing 
surface, and the parameter h denotes the metric 
coefficients, with e the angle between the 
coordinate lines x = const and z = const. For an 
orthogonal system, e = */2. The parameters K] and 
K 2 are known as the geodesic curvatures of the 
curves z = const and x = const, respectively. 
Equations (1) to (3) are subject to the follow- 
ing boundary conditions 


y ■ 

0, u = 0, 

v * 0, w * 0 

(4a) 

y * 

5, u = u e 

(x,z) , w = w e (x,z) 

(4b) 

The 

solution 

of the above equations 

also 


requires Initial conditions on two Intersecting 
planes; those In the (y,z) plane at a specified 
chordwlse location are determined from the solu- 
tions of the equations discussed In Subsection 
4.3. Those on the (x,y) plane, at a specified 
spanwlse location z » z 0 , with z 0 corresponding 
to, say, the root location, are determined from 
the solutions of the quasl-three-dlmenslonal form 
of Eqs. (1) to (3) with all derivatives with 
respect to z neglected, that Is, 


U 3u 
h^ 3x 


§7 (uh^slne) + (vh^ sine) = 0 


*■ v - K 1 u 2 cote + K 2 w 2 cosece + K^uw 


cosec e 3£ a au, 
Ph, ax + v 77 (b P 


(5) 


( 6 ) 


U 3W 

h^ + V 


3W 

ay 


K^w cote ♦ 


2 

u cosece ♦ 


21 


uw 


cote cosece ae . a /h aw. 
Ph 1 ax * u 57 < b »y) 


( 7 ) 


subject to the same boundary conditions given In 
Eq. (4). 


3.2 Interaction Law 

To account for possible flow separation, as In 
two-dimensional flows, we use the Interaction law 
of Veldman 13 where, for airfoil flows, the edge 
velocity Is expressed as the sum of an Invlscld 
velocity ug( x) and perturbation velocity 4u e (x) 
due to viscous effects, that Is, 


u e (x) = uji(x) + 4u e (x) (8) 

The perturbation velocity Is given by the Hilbert 
Integral 


“.<«> • ; i 6 h 

x a 


(9) 


In the Interaction region (x a , x b ). 


To extend this Inverse formulation to three- 
dimensional flows, It Is necessary that the two- 
dimensional Interaction formula given by Eq. (9) 
be either modified to account for the Interaction 
in the x- and z-dlrectlons or be replaced by 
another formulation which provides a relationship 
between displacement surface and external veloc- 
ity. Here we use the former approach, as des- 
cribed In Ref. 14, and first generate an Initial 
displacement surface by solving the quasl-three- 
dlmenslonal boundary-layer equations subject to 
the boundary conditions given by Eqs. (4) and (8) 
with the external velocity distribution u§( x) 
obtained from the panel method. The second step 
Involves Interaction between the Invlscld flow 
equations and the quasl-three-dlmenslonal flow 
equations. As In two-dimensional flows, the 
solutions of the boundary-layer equations are used 
to compute distributions of blowing velocity on 
the surface, and these allow the Invlscld flow 
solutions to be updated. In step three, after the 
calculation of the Initial conditions In the (y,z) 
and (x,y) planes, the fully three-dimensional 
boundary-layer equations are solved with the 
external velocity components resulting from step 
two. As before, the spanwlse velocity component 
Is assumed to correspond to Its Invlscld value. 
The viscous flow solutions are obtained by march- 
ing In the spanwlse direction at each advancing 
chordwlse location. This represents the first 
phase In an Interactive loop that Involves the 
fully three-dimensional boundary-layer equations. 
In the subsequent phases, as before, the blowing 
velocity distribution Is used to obtain Improved 
Invlscld flow solutions, so the fully three- 
dimensional boundary-layer equations are solved 
for Iced wings as for clean wings In transonic 
flow. 14 


The viscous effects in the spanwlse component 
w e are assumed to be second order, although their 
neglect Is contrary to the 1 rrotat lonal Ity condi- 
tion. However, trial calculations Involving vari- 
ations of both velocity conditions showed that 
errors were smaller than those associated with 
the convergence of the solutions. 
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3.3 Transformed Equations 


As In two-dimensional flows, we express the 
boundary-layer equations In transformed variables 
for computational purposes. At first, when the 
equations are solved for a prescribed external 
velocity distribution (standard problem), we use 
the Falkner-Skan transformation and a modified 
version of this transformation for the Inverse 
mode. In the former case, the Independent vari- 
ables are defined by 

u e 1/2 x 

x = x, z « z, n = l~) y. s = J h,dx 

vS 0 1 (10) 

For the dependent variables u, v and w, we Intro- 
duce a two-component vector potential such that 

uh^sine = |^ f wh^slne * 

»",V“ ■ - <1! • If' 

In addition, dimensionless parameters f and g are 
defined by 
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The coefficients m-j 

to are 

defined in 

Ref. 9. 

In 

terms 

of transformed variables. 

boundary 


conditions given by Eq. (4) become 

n = 0 : f =*g*f‘ *9* *0 

(16) 



The form of the transformed qua si -three-dimen- 
sional equations Is Identical to the form of Eqs. 
(13) to (15), except for small differences 
discussed In Ref. 9. 

To solve the equations in the Inverse mode, we 
define Independent variables by 

u 1/2 x 

x » x, z * z , Y - ( — ) y, S * J h dx (17) 

us o 

and relate the vector potentials t and $ to f and 
9 by 

* * (u Q vs) 1/2 h ? slnef(x,z,n) 

(18) 


* = (u Q vs) 1/2 h^slnegfx.z.n) 


and with a prime now denoting differentiation with 
respect to Y and u e and w e denoting edge velocity 
components normalized by reference velocity u 0 , 
Eqs. (1) to (3), with e’ defined by Eq. (15) and 
m‘s given In Ref. 9, are written as 

(bf-) ' ♦ ef" + m 2 [(f') 2 - (u e ) 2 ] ♦ m 5 [f‘g' - u e w e ] 
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The transformed boundary conditions for the sys- 
tem of Eqs. (19) and (20), with u e given by Eq. 
(B) and with w e corresponding to its Invlscid 
value, are 

t] * 0: f . g . f' . g' . 0 (21a) 

71 * n e : ^ 9* * ^ (21b) 


The quas 1-three-dlmensional form of the equa- 
tions, which are subject to the boundary condi- 
tions given by Eq. (21), are obtained from the 
above equations by setting 


3f 1 
3Z 


as! 

dZ 




0 and m, =0 (22) 
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To generate the initial conditions near the 
leading edge of the iced wing, we use quasl-three- 
dimenslonal boundary- layer equations expressed in 
the inverse mode given by 


(bf")”* ef" + m 2 [(f ') 2 - (u e ) 2 ] + m 5 (f , g" 

♦ V (9 ' )Z - < 5 e > 2 ] * "W fl f r- 
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The above equations can be further simplified 
If we assume that two adjacent defining sections 
of a wing are connected by straight line develop 
ment, as commonly used In the wing design. This 


j 


6 


feature simplifies the problem of shaping the 
metal for a wing surface. As a consequence, we 
can neglect the geodesic curvature of x = constant 
lines, namely K 2 , and thus set m 3 =« m 8 * 0. From 
the definitions of 1*4 and 1115 , It can be seen that 
as a result of the above assumption, these two 
terms are also small and can be neglected. We 
further assume that the local variations In cross 
sections In the spanwlse direction are small. 
Examination of the terms nvj , mg and mo for a typi- 
cal wing shows that m 2 reaches a value less than 
0.1 very close to the leading edge (x/c < 0 . 01 ) 
and mg reaches a maximum value of 0.2. However, 
their magnitudes rapidly decrease with Increasing 
x and reach a very small value at x/c < 0.1. This 
behavior allows us to neglect m 2 and mg In the 
equations and set m-j = 1 / 2 . 

3.4 Solution Procedure 

A detailed description of the solution <proced- 
ure will be reported separately. Briefly, the 
boundary-layer equations expressed In terms of 
transformed variables are solved with Keller's 
two-point finite-difference method 15 (box scheme) 
with boundary conditions expressed In Inverse form 
with the Interaction law described In Subsection 
3.2. Depending on the complexity of the flow- 
field, two forms of the box scheme are employed. 
In regions where all velocity components are pos- 
itive, the regular box scheme Is used. In regions 
of either a negative spanwlse velocity component 
or negative streamwlse velocity component, the 
zig-zag box scheme described In Ref. 15 Is used. 

3.5 Performance Degradation of an Iced Tapered 
Wing 

The Interactive boundary-layer method of Sec- 
tion 3 was used to study the performance degrada- 
tion of an Iced wing having MS-31 7 airfoil stream- 
wise sections, an aspect ratio of 3 . 43 , and a 
taper ratio of 0.4. Icing conditions were chosen 
to correspond to those In Runs 8 and 11, shown In 
Table 1 . The pressure distribution on the wing 
was computed at four locations defined by the 
midsection of each wing-section with a hundred 
panels on each defining airfoil section. The 
Ice shapes corresponding to this pressure distri- 
bution were computed with the method of Section 2 
In the middle of each wing section and were used 
to distribute Ice along the leading edge of the 
tapered wing. The computed Ice shapes for a = 
2 ° were then assumed to be the same for all angles 
of attack on the wing In the Investigation of the 
performance degradation of the wing due to Ice 
shapes corresponding to the atmospheric conditions 
given In Runs 8 and 11. At a specified angle of 
attack, with the defined Ice shapes along the 
leading-edge of the wing, calculations were per- 
formed with the method of Section 3 ; that Is, 
Invlscld flow calculations performed for an Iced 
wing were followed by the Inverse three- 
dimensional boundary-layer calculations to deter- 
mine the blowing velocity distribution to be used 
In the Incorporation of viscous effects Into the 
Invlscld method. The Invlscld flow solutions made 
use of four lifting strips, and the viscous flow 
calculations Included boundary-layer calculations 
on the wing and In the wake, the latter requiring 
velocities at off-body points In the potential 
field. This Interactive and Iterative procedure 
was repeated until the solutions converged. The 
lift coefficients were then calculated from the 


Invlscld method for each Individual strip and 
Included the contribution of Ice protruding beyond 
the wing contour and the drag coefficients from 
the boundary-layer calculations. 

Figure 6 shows the variation of the calculated 
lift coefficients as a function of angle of 
attack. Since the primary purpose of the calcu- 
lations was to demonstrate the Increase In drag 
due to Ice on a tapered wing, the angle of attack 
range was not extended to stall, which would occur 
at relatively high angles of attack for low aspect 
ratio wings. The higher lift coefficient than for 
the clean wing shown for the two Iced wings Is due 
to the normalization with the wing area of the 
clean wing In both cases. The conclusion from 
this figure Is that lift Is not affected by the 
rime Ice accretion for the angle of attack range 
considered here because the Ice shapes along the 
leading edge of the wing for runs 8 and 11 do not 
cause premature flow separation on the wing. 



Fig. 6 . Effect of leading-edge 390 and 1164 sec- 
ond rime Ice on the lift coefficient of a tapered 
wing for R = 4.6 x 10^ based on root chord. 

The calculated drag coefficients shown In Fig. 
7 represent the profile drag of the wing only and 
do not represent the total drag, since that 
requires the contribution of the Induced drag. 
The profile drag was calculated sectlonwlse from 
the Squire-Young formula based on the resultant 
velocity at the trailing edge. Comparable results 
were also obtained from the momentum deficiency 
In the far wake. Here we see considerable dif- 
ferences between the clean wing and the two Iced 
wings because the Reynolds number Is relatively 
low (Re = 4.6 x 10® for the root chord) and 



Fig. 7. Effect of leading edge 390 and 1164-second 
rime Ice on the profile drag coefficient of a 
tapered wing for R c * 4.6 x 10 5 based on root 
chord. 
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there are large regions of laminar flow on the 

clean wing. The principal contributor to the drag 
Increase for the iced wing Is the shift In tran- 
sition to near the leading edge due to roughness 
of the Iced surface. The contribution of the 
surface roughness Itself to the drag Is very small 
for Run 8 because the extent of Ice Is small and 
Its shape emulates an airfoil leading edge. The 
additional drag Increase for Run 11 results from 

the surface roughness spread over a large wetted 
area Increment, The main conclusion that can be 
drawn from these comparisons Is that drag Incre- 
ments obtained between clean and Iced airfoils In 
wind tunnels depend on transition locations on the 
clean wing. If the Run 8 case represents a wing 

with transition fixed at the leading edge and the 

clean wing case Is transition free, the observed 
drag Increments from the Run 11 case are quite 
different from each other. As a corollary, drag 
Increments obtained from wind-tunnel tests may be 
meaningless without fixing transition or knowing 
where transition occurs during the tests. 


3.6 Results for an Unswept NACA 0012 Ming 

Bragg et al. 16 have tested two wings with a 
simulated Ice shape to determine Its effect on wing 
aerodynamic characteristics. They also tested the 
same wings In clean conditions to establish the 
base case. Their measurements Include selected 
chordwlse pressure distributions, balance data on 
lift and drag coefficients, and section drag data 
by wake measurements. Since these measurements 
were conducted In a wind tunnel and our calcula- 
tions were to be done for freestream conditions, 
at first we decided to perform the calculations for 
the clean unswept wing with the Interactive method 
of Section 3. The Invlscld code used seven lifting 
strips, each with 180 chordwlse panels along the 
semi span. 

Figure 8 shows the calculated pressure distri- 
butions for a - 4 degrees together with the 
experimental results. The overall agreement 1$ 
very good. Also shown are the Integrated sectional 
lift coefficients which differ somewhat from case 
to case, but this Is expected from Integration of 
nonsmooth data. Studies are In progress to evalu- 
ate the Interactive method for the swept clean wing 
and then apply the method to both unswept and swept 
wings with Ice. 


4.0 Calculation of Ice Shapes on 
Multielement Airfoils 

To extend the method developed for analyzing 
Iced airfoils and wings to high- 1 1 ft configura- 
tions, our studies first concentrated on the calc- 
ulation of Ice shapes on the slat of a four-element 
airfoil shown In Fig. 9. Figure 10 shows the 
Invlscld pressure distribution of the clean four- 
element airfoil at a =* 0°. The Ice shapes of the 
first element corresponding to times up to two 
minutes are shown In Fig. 11 for a time step At 
of one minute. With Ice build-up on the first 
element, the computed pressure distribution 
build- up remains essentially the same except 
along the first element. Figure 12 shows the 
progression of the pressure distributions of the 
first element with time. As can be seen, the Ice 
accretions cause rapid changes In the pressure 
distribution with large leading-edge peaks. 



Fig. 8. Comparison of calculated (solid lines) and 
measured (symbols) results for the unswept clean 
wing of Ref. 16 at a = 4°: (a) y/b = 0.168, (b) 
y/b = 0.336, (c) y/b = 0.497. 
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Fig. 9. Geometry of the four-element airfoil. 
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the four-element airfoil at a = 0°. 



Mg. 12. Pressure distribution along the first 
element of the four-element airfoil with glaze Ice 
accretion at a = 0° . 

Very recently, the above method has also been 
extended to multielement airfoils. Figures 13 to 
16 show the pressure distributions and the Ice 
shapes on the first two elements of the four- 
element airfoil at a = 4° and 6°. The Ice 

shapes correspond to 2 and 5 minutes. 

Figures 17 and 18 show the results for the 
four-element airfoil at a = 0° for a two-minute 
glaze Ice computed by the multielement LEWICE code. 
Additional studies are In progress and will be 
reported later. 
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Fig. 13. Pressure distribution for the clean four- 
element airfoil at a = 4° . 
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Mg. 14. Glaze Ice shapes on the first two ele- 
ments of the four-element airfoil at a * 4°. 
The Ice shapes correspond to 2 and 5 minutes. 
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Fig. 15. Pressure distribution for the clean four- 
element airfoil at a = 6°. 



-4 0 00 4.0 60 120 

X 

Fig. 16. Glaze Ice shapes on the first two ele- 
ments of the four-element airfoil at a = 6°. 

The Ice shapes correspond to 2 and 5 minutes. 
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Fig. 17. Pressure distribution for clean and Iced 
four-element airfoil at a = 0°. 
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Fig. 18. Computed two-minute glaze Ice shapes on 
a four-element airfoil at a = 0®. 


5.0 Concluding Remarks 

Until recently, the only capability for pre- 
dicting Ice shapes on aerodynamic configurations 
was limited to single airfoils. With the methods 
described here and with the method described In 
Ref. 17 for wings, this capability now Includes 
wings and multielement airfoils. These methods, 
however, are In their Infancy and require Improve* 
ments and validation with experimental data. 
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The Interactive method for three-dimensional 
flows also provides a new capability that, except 
for the recent work of Ref. 18, did not exist for 
Iced wings. Both methods are, also In the early 
development stages and require additional work and 
validation with experimental data. 
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